clear all
close all
% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the standard limited commitment village
%        economy for both general and partial equilibrium

% This version compares two income processes

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version September 2021
% =====================================================================

global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table AGG vs vss ScaleRov

filetosave='Vilcoal_replicate_June_2018';
write=1; % set this to 1 if you want to save tables, workspace and figures
%====================================================
% I. General equilibrium
% ====================================================
% Here I use her estimates for Aurepalle village, setting heterogeneity
% parameters to average values
nyvec=[3,6]; % number of indiv income values
rho=1.0;
server=0;
if server==0
    datapath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Data\Output';
    programpath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes';
else
    datapath='/home/workgroups/testob';
end
relvarincscale=1;  % determines how to scale the difference in consumption growth variance
analytic =1; % analytic VCV used - not needed
Pun_fac=0.0;
datafile='armoments_fe0';
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);


COUNT=0;
for village=1:3
    cd(programpath)
    rho1=covyylag(village)/vary(village) % AR(1) of income
    vareps=(vary(village)).*(1-rho1.^2);
    if village==1
        N=34-1;                 % the number of total village members is N+1 under individual deviations
        beta=0.89;
        
    elseif village==2;
        N=37-1;                 % the number of total village members is N+1 under individual deviations
        beta=0.90;
    elseif village==3
        N=31-1;                 % the number of total village members is N+1 under individual deviations
        beta=0.88;
    end
    Nvill=1;
    Nincproc=1;         varnu=0;                % number of income processes to look at
    incproc=1;
    T=6200; Tinit=201;                 % number of simulation periods in total, and number of initial periods to discard
    maxcerr=0;                          % max cons meas error in log levels
    N_cerr=2;                           % number of support points for cons meas error
    maxincerr=2/3;                      % max inc meas error in multiples of inc variance (in levels)
    TT=50000; % number of periods
    vss=N;
    nyagg=5; % number of agg inc values
    Qrho=1; % number of rho parameters
    Qbeta=1;% number of beta parameters
    % set this to 1 if you want to calculate the moments for levels of
    % consumption based on log consumption
    momentclevinlog=1;
    
    for ny=nyvec % loop over income processes
        COUNT=COUNT+1;
        % 1. Income Process
        % specify a markov process without persistence for individual income
        % now we use Rouwenhorst
        [incomevalsalt, Q1] = rouwenhorst(rho1,vareps^0.5,ny);
        incomevals=reshape(exp(incomevalsalt'),length(incomevalsalt),1);
        
        %  Simulate HH income for N HH
        QQ1=Q1^1000;
        SSdist=QQ1(1,:);
        Qsum=cumsum(Q1);
        for jj=1:N
            z(jj,1)=min(find(rand(1,1)<=cumsum(SSdist)));
            y_sim(jj,1)=incomevals(z(jj,1));
        end
        for n=1:N
            for t=2:TT
                z(jj,t)=min(find(rand(1,1)<=cumsum(Q1(z(jj,t-1),:))));
                y_sim(n,t)=incomevals(z(jj,t));
            end
        end
        
        % estimate aggregate income process
        XX=[ones(1,TT-100)',log(sum(y_sim(:,100:end-1),1))'];
        YY=log(sum(y_sim(:,101:end),1))';
        beta_rov=inv(XX'*XX)*XX'*YY;
        eps_Y=(YY'-beta_rov'*XX')';
        Var_eps_Y=(eps_Y'*eps_Y)/TT;
        [S_rov_temp, P_rov] =rouwenhorst(beta_rov(2),Var_eps_Y^0.5,nyagg);
        S_rov=exp(beta_rov(1)/(1-beta_rov(2))+S_rov_temp)';
        S_rov=S_rov/ScaleRov;
        
        % ======== c. state of the economy: combinations of indiv and aggregate income ============
        S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
        P=kron(Q1,P_rov);
        cumP=cumsum(P')';
        cumQ=cumsum(Q1')';
        state=length(S);
        Y=S(:,1)+ScaleRov*S(:,2);%sum(S,2);				% possible aggregate income outcome
        % ============ transition matrix for state of the economy ============
        agent=2;
        clear caut eeewaut waut ewperf share normperf eewperf wsbnewsol wsb
        % Compute set of autarky values
        for k=1:agent
            caut(:,k)=S(:,k)*(1-Pun_fac);
        end
        waut=zeros(state,agent);
        % autarky values
        eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
        eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
        % perfect insurance values
        ewperf(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,Y./(vss+1));
        share=linspace(0.3/(1+ScaleRov),2/(1+ScaleRov),100); % constant share of individual
        for i=1:100
            ewperf(:,1,i)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,Y*share(i));
            ewperf(:,2,i)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,Y*(1-share(i))/ScaleRov);
            normperf(i)=max(max(abs((eeewaut(:,:)-ewperf(:,:,i)))));
        end
        % check if perfect insurance is feasible
        [mintemp,minind]=min(normperf);
        eewperf=ewperf(:,:,minind);
        clear ewperf mintemp minind
        if eewperf-eeewaut>0
            disp('perfect insurance')
            %  pause
        end
        
        for j=1:size(S,1)
            waut(j,1)=utility(rho,caut(j,1))+beta*P(j,:)*eeewaut(:,1);
            % autarky values for the rest of village are always the same - perfect
            % risk-sharing among the rest of village
            waut(j,2)=utility(rho,caut(j,2))+beta*P(j,:)*eeewaut(:,2);
            
        end
        
        % ====================================================
        % II. Compute constrained insurance
        % ====================================================
        
        % grid of time-varying planner weights
        lambda=[0.005,0.1,0.2,0.25,0.3,0.35:0.05:0.65,0.7,0.75,0.8,0.9,0.995]';
        m=length(lambda);
        gammagrid=[lambda vss/ScaleRov*(1-lambda)];
        gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
        
        % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
        n=length(gammagrid);
        w=zeros(n,state,agent);
        gammar=zeros(n,state,agent);
        c=zeros(n,state,agent);
        phi=zeros(n,state,agent);
        
        % consumption given aggregate resources and planner weights
        for j=1:state
            % for k=1:agent
            c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
            c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
            % end
        end
        cfirstbest=c;
        %  relative planner weights
        for j=1:state
            for k=1:agent
                gammar(:,j,k)=gammagrid(:,agent+k);
            end
        end
        
        % Compute the first best value functions without binding
        % constraints
        w=zeros(n,state,agent);
        %T first best expected utility
        for j=1:state
            w1(:,j,1)=utility(rho,cfirstbest(:,j,1))+beta/(1-beta)*[utility(rho,cfirstbest(:,:,1))]*P(j,:)';
            w1(:,j,2)=utility(rho,cfirstbest(:,j,2))+beta/(1-beta)*[utility(rho,cfirstbest(:,:,2))]*P(j,:)';
        end
        dww=1;
        wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
        while dww>0.00001
            % w(:,:,1)=utility(rho,cfirstbest(:,:,1))*inv(diag(ones(1,state),0)-beta.*P)
            %w(:,:,2)=utility(rho,cfirstbest(:,:,2))*inv(diag(ones(1,state),0)-beta.*P)
            w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
            w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
            dww=max(max(max(abs(w-wold))));
            wold=w;
        end
        wfirstbest=w;
        
        % continuation utility as function of planner weight and state of the
        % economy
        for k=1:m
            for j=1:state
                ewfirstbest(k,j,1)=P(j,:)*wfirstbest(k,:,1)';
                ewfirstbest(k,j,2)=P(j,:)*wfirstbest(k,:,2)';
            end
        end
        for k=1:m
            for j=1:state
                ewautaut(k,j,1)=P(j,:)*waut(:,1);
                ewautaut(k,j,2)=P(j,:)*waut(:,2);
                wautaut(k,j,1)=waut(j,1);
                wautaut(k,j,2)=waut(j,2);
            end
        end
        
        % ====== a. Start the loop with the enforcement constraints ======
        %%this will be the parameterized expectation
        ewsb=ewfirstbest;
        wsb=wfirstbest;
        gap=1;
        options=optimset('Display','off');
        count=0;
        eps=0.0000001;
        t=0;
        ggammagrid=fix(gammagrid*100);
        itgap=0;
        maxitgap=200;
        while gap >=0.01 && itgap<=maxitgap
            clear rr
            itgap=itgap+1;
            coalnonsust=zeros(n,state);
            count=0;
            % Begin proper loop
            count=0;
            countgrid=[ones(n,state)];
            for j=1:state
                %no permutations here
                disp(j)
                for i=1:n
                    count=count+1;
                    index=[i,j];
                    
                    % Now start updating all the functions c, z, w
                    % 1. Constraint binding for Person 1 only
                    if utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)<waut(j,1) && utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)>=waut(j,2) ;
                        bind=1; %i.e. agent 1 has the binding constraint
                        rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1));
                        if isempty(rr)==1
                            coalnonsust(i,j)=1;
                            nonconv(i,j)=nan;
                        else
                            %x0 is the guess for consumption of constrained agent
                            %1, corresponding to planner weight rr(1)
                            x0(1) = cfirstbest(rr(1),j,1);
                            %vectors of consumption and utiltiy that meet
                            %participation constraints
                            
                            [x, output]=agent1maxrho_rov_15_Mar(x0,index,ewsb,waut,bind);
                            
                            %This keeps track of nonconvergence, but
                            %allows the programme to continue - sometimes it
                            %converges in a subsequent iteration
                            if output(3) ~=1
                                %disp('Bind_ind in case j,i, did not converge')
                                % disp([j,i,output(1)-waut(j,1),output(3)])
                                nonconv(i,j)=output(3);
                                wsbnewsol(i,j,1:agent)=wsb(i,j,1:agent);
                            else
                                nonconv(i,j)=nan;
                                csol(i,j,1)=x(1);
                                gammar(i,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                csol(i,j,2)=(Y(j)-x(1))/ScaleRov;
                                wsbnewsol(i,j,1:agent)=output(1:agent);
                                phisol(i,j,1)=gammagrid(i,4)^rho/gammar(i,j,2)^rho-1;
                                phisol(i,j,2:agent)=0	;
                            end
                            %updated Lagrange multipliers using Ligon,
                            %Thomas and Worrall formulation and consumption shares instead of Pareto weights (cf LTW eq 17)
                        end
                        
                        % 2. Constraint binding for person 2 only
                    elseif utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)>=waut(j,1) && utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)<waut(j,2);
                        bind=2; %i.e. agent 2 has the binding constraint
                        rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                        if isempty(rr)==1
                            coalnonsust(i,j)=1;
                            nonconv(i,j)=nan;
                        else
                            x0(1) = cfirstbest(rr(end),j,1);
                            [x, output]=agent1maxrho_rov_15_Mar(x0,index,ewsb,waut,bind);
                            
                            if output(3)~=1
                                %disp('Bind_rov in j,i did not converge')
                                %disp([j,i,output(2)-waut(j,2),output(3)])
                                nonconv(i,j)=output(3);
                                wsbnewsol(i,j,1:agent)=wsb(i,j,1:agent);
                            else
                                nonconv(i,j)=nan;
                                csol(i,j,1)=x(1);
                                gammar(i,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                csol(i,j,2)=(Y(j)-x(1))/ScaleRov;
                                wsbnewsol(i,j,1:agent)=output(1:agent);
                                phisol(i,j,2)=gammar(i,j,2)^rho/gammagrid(i,4)^rho-1;
                                phisol(i,j,1)=0;
                            end
                        end
                        
                        % 3. Constraint binding for none
                    elseif utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)>=waut(j,1) && utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)>=waut(j,2) ;
                        csol(i,j,1:agent)=cfirstbest(i,j,1:agent);
                        gammar(i,j,1:agent)=gammagrid(i,agent+1:2*agent);
                        phisol(i,j,1:agent)=0;
                        for g=1:agent
                            wint(g)=interp1(gammagrid(:,4),ewsb(:,j,g),gammagrid(i,4));
                        end
                        wsbnewsol(i,j,1)=utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1);
                        wsbnewsol(i,j,2)=utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2);
                        nonconv(i,j)=nan;
                    end % of if loop
                end  % end of i loop                
            end  % end of j loop
            
            if sum(sum(coalnonsust))>0 || sum(sum(isnan(nonconv)))<state*n  && itgap>=3
                % break
                % disp('break')
                coalnonsustind(Qrho,Qbeta,vss)=mean(mean(coalnonsust(find(coalnonsust>0))));
                nonconvind1(Qrho,Qbeta,vss)=sum(sum(isnan(nonconv)));
                gap=100;
            else
                %T 29_Mar is the norm taken over agent 1 only?
                gap=norm(wsbnewsol(:,:,1)-wsb(:,:,1))
            end
            nonconvind=mean(nonconv(find(~isnan(nonconv))));
            if min(min(nonconv))<0
                disp('Problem with convergence in this iteration')
                %pause
                %Nonconvind(vss,Qrho,Qbeta)=mean(nonconv(find(nonconv(:,:)>0,vss)))
            end
            wsb=wsbnewsol;
            for k=1:m
                for j=1:state
                    ewsb(k,j,1)=P(j,:)*wsb(k,:,1)';
                    ewsb(k,j,2)=P(j,:)*wsb(k,:,2)';
                end
            end
        end % of iteration loop
        
        % =========================================================
        % Simulate the economy
        % =========================================================
        clear csim Yagg income indSjj incomeind gammarfun phinew
        csim=nan(vss,T);
        Yagg =nan(1,T);
        income =nan(vss+1,T);
        t=1;
        %initial state
        indSjj(:,1)=randi(ny,N+1,1);
        indSjj(:,1)=randi(ny,N+1,1);
        stream = RandStream.getGlobalStream;
        for k=1:vss+1
            incomeind(k,t)=randi(ny);
            income(k,t)=incomevals(incomeind(k,t));
            csim(k,t)=income(k,t);
        end
        Yagg(t)=sum(income(:,t));
        [Smintemp(t),Sind(t)]=min(abs(Yagg(t)/(vss+1)-S(1:nyagg,2)));
        for k=1:vss+1
            indSjj(k,t)= (incomeind(k,t)-1)*nyagg+Sind(t);
        end
        
        while t<T
            t=t+1
            for k=1:vss+1
                incomeind(k,t)=min(find(cumQ(incomeind(k,t-1),:)>=rand(1,1)));
                income(k,t)=incomevals(incomeind(k,t));
            end
            Yagg(t)=sum(income(:,t));
            [Smintemp(t),Sind(t)]=min(abs(Yagg(t)/(vss+1)-S(1:nyagg,2)));
            for k=1:vss+1
                indSjj(k,t)= (incomeind(k,t)-1)*nyagg+Sind(t);
            end
            
            for k=1:vss+1
                % calculate weights at beginning of period
                gammarfun(k,t)=(Yagg(1,t-1)-csim(k,t-1))/csim(k,t-1)/ScaleRov;
                gammarfun(k,t)=min(gammarfun(k,t),max(gammagrid(:,4)));
                % calculate the new planner weight for all individuals
                phinew(k,t)=interp1(gammagrid(:,4),phisol(:,indSjj(k,t),1),gammarfun(k,t));
            end
            
            phinew(:,t)';
            x0=csim(:,t-1)*Yagg(1,t)/Yagg(1,t-1);
            vs=vss+1;
            % solve for consumption
            [x,output]=conssolve(x0,phinew(:,t),csim(:,t-1),Yagg(1,t));
            csim(:,t)=x';
            if isnan(csim(:,t))
                stop
            end
        end
        
        
        % =========================================================
        % calculate moments of consumption and income growth
        % =========================================================
        clear incgrowth cgrowth cshare_raw dcshare_raw
        csim0=csim;income0=income;
        csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
        cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
        incerr=randn(size(csim0));   Momentsest1=nan(30,1); ind=[];
        % loop over size of cons measurement error
        for kkk=1:size(csim0,1)
            tempvar(kkk)=var(log(csim0(kkk,:)));
        end
        Tempvar=mean(tempvar);
        ic=1;
        csim=exp(log(csim0));
        logcsim=log(csim);
        logcsim0=log(csim0);
        logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
        logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
        % add inc meas error back in
        income=exp(log(income0)+varnu(incproc)^0.5*incerr);
        logincome=log(income);
        logincome0=log(income0);
        logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
        income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
        
        for ivill=1:Nvill
            aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
            logaggvillincome=log(aggvillincome);
            logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
            logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
            logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
        end
        aggincome=sum(income,1);
        
        % === discard first Tinit-1 periods
        % this chooses log or level consumption for
        % moments
        if momentclevinlog==1
            ccc=logcsim(:,Tinit:T);
            ccc_cond=logcsim_meandev(:,Tinit:T);
            ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
        else
            ccc=(csim(:,Tinit:T));
            ccc_cond=(csim_meandev(:,Tinit:T)); %NB this doesn't actually exist...
        end
        yyy=logincome(:,Tinit:T);
        yyy_cond=logincome_meandev(:,Tinit:T);
        yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
        cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
        cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
        cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
        yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
        yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
        yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
        aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
        aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
        vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
        
        % calculate moments for each individiual
        for ind=1:size(csim,1)
            dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
            dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
            dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
            dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
            
            % conditioning on income increases, and income
            % decreases
            % NB: we allocate no change to decrease!!
            rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
            rrhc=find(yy(ind,:)>0);
            rrlc=find(yy(ind,:)<=0);
            rrhc_cond=find(yy_cond(ind,:)>0);
            rrlc_cond=find(yy_cond(ind,:)<=0);
            rrhc_groupcond=find(yy_groupcond(ind,:)>0);
            rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
            
            
            %THE VARIANCES
            vardclog(ind,vss)=var(cc(ind,:));
            vardylog(ind,vss)=var(yy(ind,:));
            vardc_dypos(ind,vss)=var(cc(ind,rrhc));
            vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
            
            %conditional on village income!
            vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
            vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
            vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
            vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
            
            %%conditional on group income!
            %Step 1: create your residuals
            vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
            vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
            vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
            vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
            vardc_dypos(ind,vss)=var(cc(ind,rrhc));
            vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
            
            %%and the relative variances
            reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
            reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
            reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
            
            
            %%THE BETA COEFFICIENTS
            %1) Unconditional
            %risk-sharing whole sample
            covtemp=cov(cc(ind,:),aggyy(1,:));
            betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
            covtemp=cov(cc(ind,:),yy(ind,:));
            betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
            %risk-sharing for those with income gains versus
            %those with income losses
            YY=cc(ind,:)';
            yy_hc(ind,:)=zeros(1,size(cc,2));
            yy_hc(ind,rrhc)=yy(ind,rrhc);
            oneone=ones(1,size(cc,2));
            ones_hc(ind,:)=zeros(1,size(cc,2));
            ones_hc(ind,rrhc)=oneone(rrhc);
            XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
            beta_yhigh(:,ind,vss)=regress(YY,XX);
            betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
            betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
            
            %2) conditional on villag income
            covtemp=[];
            X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
            YY=cc_cond(ind,:)';
            beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
            
            X=[ones(size(aggyy(1,:),2),1),aggyy'];%
            YY=cc(ind,:)';
            beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
            
            %risk-sharing for those with income gains versus
            %those with income losses
            YY_cond=cc_cond(ind,:)';
            yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
            yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
            oneone=ones(1,size(cc_cond,2));
            ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
            ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
            XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
            beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);
            
            %3) conditional on group income
            covtemp=[];
            X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
            YY=cc_groupcond(ind,:)';
            beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
            X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
            YY=cc(ind,:)';
            beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
            
            %risk-sharing for those with income gains versus
            %those with income losses
            YY_groupcond=cc_groupcond(ind,:)';
            yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
            yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
            oneone=ones(1,size(cc_groupcond,2));
            ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
            ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
            XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
            beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
        end
        var_dyagglog(1,vss)=var(aggyy);
        % vector of moments
        % col 1-5 are parameters
        Momentsest1(1:5,ic)=[vss,beta,rho,rho,varnu]';
        % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
        % (8:9), conditonal on group income (10:11)
        Momentsest1(6:11,ic)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
            mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
            mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
        % 12-14: beta, beta high, beta low unconditional
        Momentsest1(12:14,ic)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
        % 13_16: beta, beta high, beta low conditional on
        % village income
        Momentsest1(15:17,ic)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
        % 18:20: beta, beta high, beta low conditional on
        % group income
        Momentsest1(18:20,ic)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
        %21:22 are variances of dclog, dylog
        Momentsest1(21:22,ic)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
        %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
        %unconditinoal 23:24, conditional on village income 25:26,
        %conditional on group income 27:28
        Momentsest1(23:24,ic)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
        Momentsest1(25:26,ic)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
        Momentsest1(27:28,ic)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];
        Momentsest1(30,ic)=[mean( vardylog_cond(1:size(csim,1),vss))];
        %%and the critical gap
        Vardylog(village)=Momentsest1(22,1);
        Vardylog_cond(village)=Momentsest1(30,1);
        
        Momentstable(1:3,COUNT)=Momentsest1(1:3,1);
        if relvarincscale==1
            Momentstable(4:7,COUNT,incproc)=[Momentsest1([8 15],1) ; Momentsest1(9,1)/Vardylog_cond(village) ;(Momentsest1([16],1)-Momentsest1([17],1))];
        elseif relvarincscale==2
            Momentstable(4:7,COUNT)=[Momentsest1([8 15],1) ; Momentsest1(9,1)/Momentsest1(8,1)/Vardylog_cond(village) ;(Momentsest1([16],1)-Momentsest1([17],1))];
        elseif relvarincscale==0
            Momentstable(4:7,COUNT)=[Momentsest1([8 15],1) ; Momentsest1(9,1) ;(Momentsest1([16],1)-Momentsest1([17],1))];
        end
        Momentstable(1,COUNT)=Momentstable(1,COUNT)+1;
        Momentsrepl(:,find(ny==nyvec))=Momentsest1(:,1);
    end
    
    cd(datapath)
     datafile=['icrisatmoments0'];
    eval(['load ' datafile '.out'])
    eval(['datamoments=' datafile ';']);
    eval(['clear ' datafile])
    
    Nom=4;
    bootstrapvar=datamoments(1:Nom,4+(village-1)*Nom+1:4+(village)*Nom);
    analyticvar=datamoments(1:Nom,16+(village-1)*Nom+1:16+(village)*Nom);
    analyticmoments=datamoments(1:Nom,1:3);
    analyticrawmoments=datamoments(1:Nom,28+1:28+3);
    bootstraprawvar=datamoments(1:Nom,28+4+(village-1)*Nom+1:28+4+(village)*Nom);
    analyticrawvar=datamoments(1:Nom,28+16+(village-1)*Nom+1:28+16+(village)*Nom);
    
    VCV2data=bootstrapvar;
    VCV2rawdata=bootstraprawvar;
    
    if analytic==1
        VCV2data(2,2)=analyticvar(2,2);
        VCV2data(4,4)=analyticvar(4,4);
        VCV2rawdata(2,2)=analyticrawvar(2,2);
        VCV2rawdata(4,4)=analyticrawvar(4,4);
    end
    Momvec2data=analyticmoments(:,village)';
    Momvec2rawdata=analyticrawmoments(:,village)';
    Momentsdatatable(:,village)=[nan;nan;nan;Momvec2data'];
end

rowLabels = {'n','$\delta$','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}}{Var_{dc|dy\leq 0}}$', '$\frac{\beta{dcdy>0}}{\beta{dcdy\leq 0}}$'};
columnLabels = {'Data','ID, ny=3', 'ID, ny=6','Data','ID, ny=3', 'ID, ny=6','Data','ID, ny=3', 'ID, ny=6'};
if write==1
    matrix2latex_newnew([Momentsdatatable(:,1) Momentstable(:,1:2) Momentsdatatable(:,2) Momentstable(:,3:4)  Momentsdatatable(:,3) Momentstable(:,5:6) ],['C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\Documents\Latex\momentsreplic_13_July_21.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', 'small');
end

save(filetosave)
